#------------------------------------------------------------------
# File:         wake_reynolds_check.py

# Purpose:      Based on Jensen's simple wake model, this
#               code checks whether the Reynolds number independance 
#               assumption can be applied in a CFD wake 
#               calculation for a given WT Thrust coefficient;
#               Power eficiency with respect to upstream
#               wind speed is also investigated.

# Input file:   wind turbine generator file (ASCII format)
#               1st column: wind speed
#               2nd column: Power 
#               3rd column Thrust 

# User input    - Rotor diameter in [m]
#               - Location: (offshore/onshore)

# Usage:        python wake_reynolds_check.py wtfile       

# Author:       Ben Martinez
# Created:      26/09/12
# Company:      Vattenfall A/B
#-----------------------------------------------------------------

import sys
import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate

fu = np.genfromtxt(sys.argv[1], dtype = None)

U_turb = [i[0] for i in fu]
Power  = [i[1] for i in fu]
Ct = [i[2] for i in fu]

D = float(input("\nEnter the Wind Turbine diameter in [m]: "))
loc = raw_input("\nOffshore or onshore site? (Off/On): ")

if loc == 'Off':
    k = 0.05
else:
    k = 0.075

for pp, nn in enumerate(Power):
    
    if nn == 0.0:
       Power[pp] = 0.00001
                 
X = [3*D,4*D,5*D,6*D,7*D,8*D,10*D]
X_str = ['X = 3*D','X = 4*D','X = 5*D','X = 6*D','X = 7*D','X = 8*D','X = 10*D']

U_turb_f = list(np.arange(0,U_turb[0])) + U_turb + \
list(np.arange(U_turb[-1]+1,U_turb[-1]+3))

Power_f = list(np.zeros(U_turb[0])) + Power + \
list(np.zeros(2))

f = interpolate.interp1d(U_turb_f,Power_f)
T_name = sys.argv[1].split('.')
T_name_n = '.'.join(T_name[0:(len(T_name)-1)])

for m, x in enumerate(X):
    
    Urat = []
    V = []

    for j in range(0,len(Ct)):

        calc = 1. - ((1. - (1.- Ct[j])**(1./2.))*((D/(D + 2.*k*x))**2.))
        Urat.append(calc) 
        V.append(calc*U_turb[j])
        
    
    Power_V = f(V)   
    Power_rat = Power_V[2:]/Power[2:]
    plt.figure(1)   
    plt.plot(U_turb,Urat,label=X_str[m])
    plt.legend(loc = 'best')
    plt.grid()
    plt.xlabel('upstream U')
    plt.ylabel('wake U / upstream U')
    plt.title(T_name_n + ' with ' + str(D) + 'm Rotor Diameter' )
    plt.hold(True)
            
    plt.figure(2)        
    plt.plot(U_turb[2:],Power_rat, label = X_str[m])
    plt.legend(loc = 'best')
    plt.grid()
    plt.xlabel('upstream U')
    plt.ylabel('Relative Power')
    plt.title(T_name_n + ' with ' + str(D) + 'm Rotor Diameter')
    plt.hold(True)

    
plt.show()




    
             
        
    
    
    







        


